
#include <iostream>
#include <fstream>
#include <iomanip>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

#define GSL_MONTE_FN_EVAL(F, x) (*((F)->f))(x, (F)->dim, (F)->params)

//#define PASS cout<<"##### PASS"<<__LINE__<<" #####"<<endl

using namespace std;

struct input
{
    int Nt;
    double dt;
    double x;
};

double intgrnd(double *x, size_t dim, void *prmtrs)
{
    struct input *input_p = (struct input *)prmtrs;

    int Nt = input_p->Nt;
    double dt = input_p->dt;
    double x_val = input_p->x;

    double sum = 0.0;
    double N_prpgtr = pow(1.0 / (2 * M_PI * dt), 0.5 * (Nt + 1));
    for (int j = 0; j < Nt - 1; j++)
    {
        sum = sum + 0.5 * (x[j + 1] - x[j]) * (x[j + 1] - x[j]) / dt + dt * 0.5 * x[j] * x[j];
    }
    sum = sum + 0.5 * (x[0] - x_val) * (x[0] - x_val) / dt + dt * (0.5 * x_val * x_val) + 0.5 * (x[Nt - 1] - x_val) * (x[Nt - 1] - x_val) / dt + dt * 0.5 * x[Nt - 1] * x[Nt - 1];

    return N_prpgtr * pow(M_E, -1 * sum);

    /****** fill the correct integrand here ******/
}

int main(int argc, char **argv)
{
    struct input inputs;
    inputs.Nt = atoi(argv[1]);
    inputs.dt = atof(argv[2]);
    double x_max = atof(argv[3]);
    double Nx = atof(argv[4]);
    size_t calls = atoi(argv[5]);

    int Nt = inputs.Nt;
    double xl[Nt] = {0.0};
    double xu[Nt] = {0.0};

    for (int i = 0; i < Nt; i++)
    {
        xl[i] = -6.0;
        xu[i] = 6.0;
    }

    const gsl_rng_type *T;
    gsl_rng *r;

    gsl_monte_function G;
    G.f = &intgrnd;
    G.dim = Nt;

    double res = 0.0;
    double err = 0.0;

    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(Nt);
    ofstream dat_file;
    dat_file.open("CPU_hrmnc_oscltr.dat", ios::app);
    for (int i = 0; i < Nx; i++)
    {
        inputs.x = i * x_max / (Nx - 1);

        G.params = &inputs;
        gsl_monte_vegas_integrate(&G, xl, xu, Nt, 10000, r, s, &res, &err);

        do
        {
            gsl_monte_vegas_integrate(&G, xl, xu, Nt, calls, r, s,
                                      &res, &err);
        } while (fabs(gsl_monte_vegas_chisq(s) - 1.0) > 0.25);

        printf("x = %.6e, \n\t%.6e +- %.6e, chisq = %.3f\n", inputs.x, res, err, gsl_monte_vegas_chisq(s));

        dat_file << setprecision(10) << inputs.x << "," << res << "," << err << "," << gsl_monte_vegas_chisq(s) << endl;
    }
    dat_file.close();
    gsl_monte_vegas_free(s);

    gsl_rng_free(r);

    return 0;
}